Accelerating material property prediction using generically complete isometry invariants

Periodic material or crystal property prediction using machine learning has grown popular in recent years as it provides a computationally efficient replacement for classical simulation methods. A crucial first step for any of these algorithms is the representation used for a periodic crystal. While similar objects like molecules and proteins have a finite number of atoms and their representation can be built based upon a finite point cloud interpretation, periodic crystals are unbounded in size, making their representation more challenging. In the present work, we adapt the Pointwise Distance Distribution (PDD), a continuous and generically complete isometry invariant for periodic point sets, as a representation for our learning algorithm. The PDD distinguished all (more than 660 thousand) periodic crystals in the Cambridge Structural Database as purely periodic sets of points without atomic types. We develop a transformer model with a modified self-attention mechanism that combines PDD with compositional information via a spatial encoding method. This model is tested on the crystals of the Materials Project and Jarvis-DFT databases and shown to produce accuracy on par with state-of-the-art methods while being several times faster in both training and prediction time.


Introduction
A solid crystalline material is made up of a periodically repeated unit cell containing a motif of atoms (ions or molecules).Crystals can distinguish themselves by atomic types (chemical elements and possibly charges of ions) and by the geometry of atomic centers.Both of these aspects can determine the various properties of a crystal.Knowledge of these properties is pertinent for determining whether a crystal can be experimentally synthesized or is useful for a particular application.
Determination of property values can be done using ab initio calculations with techniques like density functional theory (DFT) [1].These techniques are often computationally expensive [2].Further, they require extensive domain knowledge to be applied correctly, making them inaccessible.Recently, machine learning has become very popular as a substitute and has experienced success in decreasing computational costs while producing accurate predictions.
Any learning algorithm requires an input representation that adequately describes the object of interest.Objects similar to crystals, like molecules, are often treated as finite point clouds.This makes their representation more easily constructible than a representation for crystals, which are not bounded in size.
Figure 1: Classification of geometric descriptors for periodic crystals based on the properties possessed.Cell parameters consist of the unit cell lengths and angles, but this is ambiguous as there are an infinite number of unit cells.Space group is a label defined by the symmetry relations the crystal exhibits, but is sensitive to atomic perturbations.Equivariant GNNs cannot be used to distinguish periodic structures.PDF needs additional smoothing to retain continuity, introducing more parameters.The PDD is invariant, generically complete, and continuous under the EMD.
While a crystal can be described in several ways, descriptors that are easily human-interpretable, such as unit cell parameters or atomic coordinates are not useful for machine learning algorithms.Atomic coordinates do not retain invariance under rigid motion.Unit cell based descriptors are also ambiguous as there are infinitely many valid unit cells for a single structure.Such ambiguities can result in different model outputs for the same structure.Techniques such as data augmentation [3] and parameter sharing [4] can mitigate these effects but still do not guarantee the aforementioned consistency.
The structure-property relationship [5] dictates that changes in the structure of a material result in changes in its properties.Distinction between crystals then allows for distinction between their respective property values.A machine learning algorithm (for a regression task) is a map from a crystal representation to value.If a representation cannot distinguish periodic crystals then two different crystals can incorrectly be perceived to be the same and so will the output property values.Similarly, if the same crystal can be represented in different ways, consistent mapping cannot be guaranteed.
The fundamental model of a crystal is a periodic set of points at all atomic centers (even without atomic types), see details in Definition 2.1.Since crystal structures are determined in a rigid form, their strongest practical equivalence is rigid motion, which is a composition of translations and rotations in R n .We consider a slightly weaker isometry, which is a composition of rigid motion and mirror reflections.Two periodic point sets S and Q are isometric if they are related by an isometry f : R n → R n , so f (S) = Q.This work exclusively applies these concepts to dimension n = 3.
An isometry invariant I is a property (descriptor or function) that is preserved under any isometry.The values of I should be simpler than the initial periodic set, for example, a scalar, vector, or matrix.Not all invariants can be considered equally useful, however.Space groups, for example, reflect the symmetry of a given material but tiny perturbations of atomic coordinates can change the material's space group.Hence the important question is to quantify how different crystals are. Figure 1 illustrates the relationship between geometric descriptors based on their properties.A practically useful invariant should satisfy the following conditions introduced by [6] Problem 1. Find a function I on all periodic point sets in R n subject to the following conditions: 1a.Invariance: If two periodic point sets, S and Q are isometric, then I(S) = I(Q).
1b. Generic Completeness: If I(S) = I(Q), then the two periodic sets are isometric.1c.Computability: the invariant I, the metric d, and the reconstruction of S can be computed in polynomial time with respect to the size of the motif of any periodic point set S.
1d. Reconstructability: any periodic set S can be fully reconstructed from its invariant I(S).1f.Lipschitz continuity: If a periodic set Q is obtained by shifting points within periodic set S by at most ϵ, then the distance between the two periodic sets can be bound according to some distance function d such that d(I(S), I(Q)) ≤ Cϵ for some fixed constant C.
In the present work, an isometry invariant called the Pointwise Distance Distribution (PDD), defined formally in Definition 2.2, which has properties 1a and 1c (see Theorem 5.1 of [6]), and satisfies 1b and 1d for any periodic set in the general position [6] along with sufficiently large k, and inclusion of the lattice (see theorem 4.4 of [6]).Conditions 1e and 1f (see Theorem 4.3 of [6]) are satisfied under the Earth Mover's Distance (EMD) between PDDs.
The contribution of this work is a Transformer model [7] which utilizes the PDD to make predictions on the properties of materials in a highly efficient manner compared to state-of-the-art models.In doing this, the gap between unambiguous crystal descriptors and machine learning models is bridged.Use of such a representation produces results on par or better than graph-based models, despite the additional structuring of data that comes with edges and edge embeddings.This model is faster in both prediction and training speed compared to two state-of-the-art models.
To prove the method's robustness, the model is applied to the crystals of the Materials Project [8] and Jarvis-DFT [9].Further experimentation on material classification, hyper-parameter sensitivity testing, and prediction of crystal properties without compositional information is included in Supplemental Materials Sections E and D.

Related Work
Early works in crystal property prediction used more classical statistical methods like kernel regression [10] before eventually moving towards deep learning [11].More recent works have shifted to Graph Neural Networks (GNN) [12,13,14,15,16,17,18,19,20,21] due to their ability to make use of structured data.Several of these focus on predicting the properties of the crystals contained within the Materials Project [8] using a multigraph representation where vertices represent atoms and edges are embedded with the pairwise distances to an atom's nearest neighbors.Some state-of-the-art models use line graphs to incorporate more geometric information like angles and dihedrals [13,22].The derived line graphs can contain significantly more vertices and edges, incurring a higher computational cost.[23] take a physics principled approach and substitute the interatomic distances for interatomic potentials and capture a crystal's periodicity using the infinite sum of these potentials.
While effective in modeling crystal structures, graphs are discontinuous under perturbations [24].Small movements in the atomic positioning can cause significant changes to the graph's topology.Some graphs are not unit cell invariant.
Due to an infinite number of possible unit cells, the graph is then reliant on the data or the cell reduction technique used.

SE(3)
-equivariant models such as Tensor Flow Networks [25], SE(3)-Transformers [26], and SE(3)-GNNs [27] impose constraints on the set of learnable functions of the network such that the output is equivariant with respect to the input points.While effective on finite point clouds, they offer no promise of completeness with periodic point sets which is necessary for distinguishing between structures.These architectures can also be beneficial when predicting properties that are equivariant with respect to rigid motion, but the crystal properties examined here are invariant to such symmetries, and thus invariance of the model output through either the input or model architecture is required.
In addition to the properties mentioned earlier, the invariant needs to be able to be adapted for a learning algorithm.Further, it needs a way to incorporate compositional information as invariants typically only consider structure.Some invariants have been adapted for use in machine learning algorithms such as symmetry functions [28,29] and Voronoi cells [30].Both of these, however, still lack continuity.The Partial Radial Distribution Function is invariant and continuous but is not complete for homometric crystals.Smooth Overlapped Atomic Positions [31] has been incorporated into models for property prediction, but is an invariant for atomic environments, not entire structures.Coulomb ma-trices use electrostatic interactions between atoms instead of Euclidean distances and are invariant and complete for molecules but have not been proven to retain these properties for the periodic case [32].Average Minimum Distance (AMD) [33] is invariant and continuous and has been used to predict lattice energies via Gaussian Process Regression [34], but is incomplete and does not currently have a way to incorporate compositional information.

Methods
A periodic crystal can be represented as a periodic point set [35] with points located at the atomic centers of the structure.They do not differentiate between the types of atoms and instead treat every point as unlabeled.A periodic point set (periodic set) can defined like so: Definition 2.1 (Periodic Point Set).For a set of n basis vectors v 1 . . .v n ∈ R n , the lattice L is formed by the integer linear combinations of these basis vectors The PDD of a periodic set is the m × (k + 1) matrix where m is the number of atoms in the motif M and k is a positive integer indicating the number of nearest neighbors to use.Each row corresponds to a point in the motif and the entries within the row consist of the Euclidean distance to each of this point's k-nearest neighbors within the entire periodic set S. The first entry of the row is assigned to be a weight equal to 1 m (the distances follow).Once the matrix is formed, rows that are the same are collapsed into a single row and their respective weights are added.Due to very small differences between rows caused by floating point arithmetic or atomic perturbations, it is common to use a tolerance, henceforth called the collapse tolerance, that allows rows with small non-zero differences (e.g. with respect to L ∞ distance) to be treated as the same.By collapsing rows in the PDD, the resulting matrix representation is always the same for a given crystal, regardless of the unit cell.Formally, Definition 2.2 (Pointwise Distance Distribution).For a periodic set S = L + M with a set of motif points M = {p 1 , . . ., p m } within a unit cell U of lattice L, the uncollapsed PDD matrix for a parameter k ∈ N + is a m × (k + 1) matrix where the i th row consists of the row weight w i = 1 m followed by the euclidean distances d 1 . . .d k from the point p i to its k-nearest neighbors such that d 1 ≤ d 2 . . .≤ d k .If a group of rows is found to be identical (or close enough using a valid distance measure within some tolerance) then the matrix rows are collapsed and the weights of the involved rows are summed.The resulting matrix will then have less than m rows.
This matrix is referred to as PDD(S; k) for a periodic set S and positive integer k.

Periodic Set Transformer
In our model, rather than being considered a matrix of values, the PDD will be considered a set of grouped atoms.A single group of atoms corresponds to the k-nearest neighbor distances in a given row within the PDD matrix.Each member of the set will carry the weight provided by the row in the PDD.Any set A can trivially be turned into a weighted set by weighing each element by 1  |A| .When the PDD is not collapsed, then there can be more than a single occurrence of any given element, making the uncollapsed PDD a multiset.Now, let A be a multiset of the form A = {a (j) i : i ∈ [1, . . ., n], j ∈ [1, . . ., n i ]} where n i is the multiplicity of element a i and a (j) i is the j th occurrence of element a i .This multiset can be turned into a weighted set by assigning each element a i with the weight ni n .We can recover the influence of multiplicity by the use of weights in our model.When a periodic crystal has its unit cell modified, the proportion of each atom is expanded or reduced.The use of weights captures this behavior in the form of a concentration or frequency.
We use an attention mechanism to find the interactions between members of the set.The rows of the PDD contain the pairwise distance information, but they do not indicate which atoms these distances correspond to.Application of the attention mechanism can help the model learn these interactions.
Let R ∈ R m×k be the PDD matrix containing m rows without the associated weight column.Let w ∈ R m×1 be the column vector containing the weights from the PDD matrix.The initial embedding is X (0) = RW d where W d ∈ R k×d is the initial trainable weight matrix.The embedding is updated according to: where i W K and V = X (0) W V and d is the embedding dimension of the weight matrices for the query, key, and value W Q , W K , and W V respectively as described in [7].The function σ is the softmax function with the PDD weights integrated into it; σ is applied to each row z of the input matrix, and i and j are used to index entries in z and w.The i th entry of the output vector is defined by: The result is passed through a single-layer perceptron SLP .The layer normalization order described by [36] is used for increased stability during training.Equation (2) describes the case for a single attention head.When multiple attention heads are used, the PDD weights are applied to each individually and the result is concatenated before being passed to the SLP like so: where h is the number of attention heads, ⊕ is the concatenation operator and Q i , K i and V i are the query, key and value for the i th head.This process is repeated l times; this determines the depth of the model.The embeddings are finally pooled into a single vector by reincorporating the PDD weights into a weighted sum of the row vectors x i of the final embedding X (l) .
This final embedding can be passed to a perceptron layer to predict the property value.
This version of self-attention can be applied to a weighted set or distribution.The weights are applied in such a way that the output of the PST is invariant to an arbitrary splitting of rows within the PDD.We provide a formal proof of this in supplementary material Section A. An overview of the Periodic Set Transformer (PST) architecture with PDD encoding (described in the next section) can be seen in Figure 2.

PDD Encoding
While structure is a powerful indicator of a crystal's properties, there may be datasets in which it is not the primary differentiator of a set of crystals.In such cases, the composition of the atoms contained within the material has a heavy influence.The previously described transformer does a good job of utilizing the structural information within the PDD but does not provide an obvious way to include atomic composition.
Transformers for natural language processing tasks use positional encoding to allow the model to distinguish the position of words within a given sentence [37].A recent transformer model, Uni-Mol [38], which performed property prediction for molecules (among other tasks), used 3D spatial encoding first proposed by [39] to give the model an understanding of each atom's position in space, relative to one another.This encoding is done at the pair level, using the Euclidean distance between atoms and a pair-type aware Gaussian kernel [40].A transformer model for finite 3D points clouds is provided by [41] via vector attention.The case for crystals is more difficult because they are not bounded in size and can exhibit many symmetries.Fortunately, by using the rows of the PDD we can distinguish each atom with structural information.We refer to this as PDD encoding.
When rows are grouped together, they are done so by having the same k-nearest neighbor distances.Though rare, it is possible for rows corresponding to different atom types to be collapsed.If this occurs, the selection of either atom type will result in information loss.To prevent this, we add the condition that the groups must be formed on the basis of having the same k-nearest neighbor distances and the same atomic species.In this case, the periodic point set has points that are labeled according to atomic type.
For a periodic set S, let P DD(S; k) be the resulting PDD matrix with parameter k.Let R be P DD(S; k) without the initial weight column and T be the matrix whose rows correspond to the vector of atomic properties used to describe the type of atom associated with each row of P DD(S; k).The initial set of embeddings for the attention mechanism is defined as X (0) = RW s + T W c where W s and W c are initial embedding weights.By starting with a linear embedding, the PDD row can be transformed to match the dimension of composition embedding.The parameter k used can then be changed as needed to include distance information from further neighbors.
3 Results and Discussion

Prediction of Materials Project Properties
The model will be applied to the data within the Materials Project.To make fair comparisons to other models we report the performance according to Matbench [42], which contains data for various crystal properties.The error rates are reported using five-fold cross-validation with standardized training and testing sets for each fold.Further, tuning is done according to the models' authors and thus our model can be compared to others more fairly.
The crystals in the Materials Project are highly diverse in composition.For all predictions, we include the composition of the crystal with PDD encoding.To incorporate this compositional information the mat2vec atomic embeddings supplied by [43] are used.The embeddings have empirically been found to produce better performance than the onehot encoded method used by CGCNN [44].They also have the added convenience of not missing any atomic property information for certain elements.
In Table 1 we report the average mean-absolute-error (MAE) across the five test sets.We include the reported accuracies of other models to allow for comparison.The selection of models aims to present a high diversity in approaches while also coming from relatively recent publications.CrabNet [44] is the only other Transformer model listed on Matbench.This model, in terms of architecture, is the most similar to the PST.Additionally, the atomic embeddings used to describe each chemical element are the same as those used in our model.The majority of models used in crystal property prediction use GNNs.Matbench features several of these, but the model that provides the best results on several properties is coGN.coGN is a GNN that includes angular and dihedral information through the use of line graphs.As such, the amount of information used is significantly more than the PST, which uses the distribution of pairwise distances.Crystal Twins (CT) [45] is a model based on the convolutional layer developed by CGCNN [12] (also used by several other models [16,46,18]) that uses self-supervised learning to create embeddings based on maximizing the similarity between augmented instances of a crystal.
The PST and CrabNet share similarities in their construction.Both use mat2vec atomic embeddings and utilize a Transformer architecture with self-attention.While CrabNet uses fractional encoding to embed the multiplicity of each element type, we opt for the PDD-weighted attention mechanism and pooling described by Equations ( 2) and ( 4).The PST also uses PDD encoding to add structural information.This could also be done for CrabNet, but the combination of both fractional and PDD encoding is not guaranteed to aid in performance and simple summation of The performance disparity between the PST and coGN on formation and band gap energy can be difficult to discern.
GNNs allow embeddings at the vertex and edge level.These embeddings can not only carry different information but allow for simultaneous updates to each of these embeddings, adding richness to the learned representation.CoGN takes this further and updates the original edges with message-passing from the derived line graph, allowing the inclusion of angular information.The edges of the line graph are further updated by its line graph, incorporating dihedral angles.These updates allow the model to learn a better representation of a crystal in latent space, which can be necessary for these larger datasets.This points to a current limitation for the PST.By using PDD encoding, we effectively limit opportunities for such updates to a single embedding representing both the atom's properties and its structural behavior.
In Table 2 the time taken to train and test the PST and coGN are listed.Our model was trained for only 250 epochs while coGN is trained for 800 epochs.While this accounts for a large portion of the disparity, the time taken to train each model per epoch is also faster for our model in all properties except the refractive index.A batch size of 32 is used for datasets containing less than 5, 000 samples.For exfoliation energy and phonon peak, PST takes 68.1% and 92.6% of the time of coGN per epoch.The larger batch size of coGN allows it to have greater GPU utilization and thus, better training efficiency.For properties with greater than 5, 000 samples, the batch size for both models is the same.For these properties, the training time per epoch is between 65 − 70% of coGN.
The prediction time disparity is more pronounced.For band gap and formation energy, the PST makes predictions approximately five times faster than coGN.Using nested line graphs introduces significant computational cost, but coGN is able to shrink the size of the graph by using the atoms in the asymmetric unit.The unit cell based approach was initially proposed by CGCNN [12] and used in several follow-up works [18,16,13,45], but the unit cell is inherently ambiguous and unnecessarily large in terms of the number of atoms needed to fully describe a crystal's symmetry.The PDD at a collapse tolerance of exactly zero will have a number of rows less than or equal to the number of atoms in the asymmetric unit.At higher tolerances, this number will be reduced until reaching the number of unique chemical elements in the crystal.

Ablation Study
In Table 3 we list the results for each "component" within the Periodic Set Transformer.In the row indicating PDD as the component, we train and test the model using only the structural information within the PDD.In the "Composition" component we pass the atomic encoding for the elements in the crystal and their concentration in the form of the PDD weights to the model without the PDD encoding.The model is run using k = 15 and a collapse tolerance of exactly zero.
By separating out each component of the model, we can interpret the importance of each to a particular property.
Properties that experience a more significant decrease in performance when the PDD encoding is not used, can be ascribed to be more dependent on structural information.In all cases, the combination of both the composition and PDD encoding results in significantly lower error rates.We can conclude that this encoding method is effective in combining the structural and compositional information of a crystal structure.In Table 4 the effect of including the PDD weight in the attention mechanism described in Equation ( 2) and in the pooling layer described by Equation ( 4) is listed.A collapse tolerance of zero is used to remove any regularization effect (further described in Section E.3).
The exclusion of weights from both the attention mechanism and pooling decreases accuracy significantly.Doing this removes all indications of multiplicity making discernment of crystals more difficult.The inclusion of the weights in the pooling layer is more impactful than when applied in the attention mechanism.The use of the weights in the pooling layer alone allows the model to perform better when the number of samples in the dataset is low.Datasets with fewer samples likely have less diversity amongst their crystals, making the need for recognizing the multiplicity of atoms less necessary.

Prediction of Jarvis-DFT Properties
The Jarvis-DFT dataset [9] is a commonly used set of materials with VASP [47] calculated properties.The list of properties computed for the materials within the dataset is more extensive than that of the Materials Project.Its inclusion provides further evidence of the robustness of the model on an even wider variety of crystal properties.
The prediction MAE produced by PST and Matformer for 12 different properties from the dataset are included in Table 5.For Matformer, we retrain the model to ensure the training and testing sets are the same.We use the default parameters for the model defined by the authors' codebase.We make one alteration to the training procedure; the number of epochs trained is reduced to 250.The number of epochs is the same as for our model.The PST outperforms Matformer in nine of the twelve properties tested.In particular, properties for which data is sparse yield results that favor the PST significantly (i.e.exfoliation energy, e ij and d ij ).Jarvis-DFT has two band gap values that are computed for its crystals, one which uses the optimized Becke88 functional (OPT) [48] and the other uses the Tran-Blaha modified Becke Johnson potential (MBJ) [49].The latter is more accurate (when compared to experimentally observed values) but also more computationally expensive.For this reason, there are significantly fewer computed values in the database.Interestingly, the PST produces a smaller error for the more accurate band gap values compared to Matformer, but a larger error for the less accurate OPT calculated values.A possible reason for this is the smaller sample size for which the PST has shown to be more effective.The disparity in performance for formation and total energy can be attributed to Matformer's architecture which uses a GNN that updated both node and edge embeddings.This additional level of expression is helpful particularly when the size of the data grows larger, though it does come with added computational cost.
Matformer has been shown to produce even better results than the previous state-of-the-art model ALIGNN [13] while taking roughly a third of the time to do both training and prediction.In Table 6, the training and prediction time for each of the properties in the Jarvis-DFT dataset is reported for the PST and Matformer.For the training time, the validation and pre-processing times are not included.The prediction time listed is the number of seconds taken to make predictions on the test set.In the closest training time comparison, the PST is still more than six times faster than Matformer.The training times for all properties fall between six and twelve times faster for the PST compared to Matformer.The performance increase can be attributed to several factors.Primarily, Matformer relies on a line graph (similar to coGN [22]) in order to update edge embeddings.While this increases the information used and leads to richer learned embeddings, the size of line graphs is considerably larger than the graph they are derived from.This, in turn, incurs a higher computational cost.
The difference in prediction times is more significant.Exfoliation energy is predicted over fifty times faster using the PST than with Matformer.This is the closest the two models perform to each other.Notably, exfoliation energy also has the fewest samples.For the bulk of the other properties, the speedup ranges between eighty and ninety times faster for the PST.

Conclusion
The PDD is a representation for periodic crystals which is invariant to rigid motion and independent of unit cell.By using weights and creating a distribution, the PDD is able to represent an infinitely spanning object by its finite forms of behavior.Further, by collapsing rows in the PDD, the resulting representation can also be much smaller in comparison to the number of atoms within the unit cell, even when the cell is reduced.
The model is applied to the crystals of the Materials Project and Jarvis-DFT on a variety of material properties.Despite using less information in the model than more commonly employed graph-based models, the PST is able to produce results on par or even exceeding that of models like coGN and Matformer while taking significantly less time to train and make predictions.

Accelerating Material Property Prediction using Generically Complete Isometry Invariants: Supplemental Material
Jonathan Balasingham, Viktor Zamaraev, and Vitaliy Kurlin

A Equivalence of Distributions in the PST
If a PDD is arbitrarily expanded or collapsed, the PST should produce the same results as these PDDs are considered equivalent.The same can be said of any input distribution.This is proven here: Lemma A.1.Let A and B be weighted multisets each containing elements from the set S = {x 1 , . . .x n }.Each element x i ∈ R 1×n occurs with multiplicity m Proof.To prove that the output of the PST is equivalent for A and B, it is sufficient to prove that the output of σ (defined in Equation ( 2)) and the pooling layer (defined in Equation ( 4)) are the same for A and B. Let , and v i = x i W V be the query, key, and value vectors produced by the weight matrices W Q , W K , W V ∈ R n×d .First, note the pre-softmax attention weight from x i to x j can expressed as: and is independent of both weight and multiplicity.Further, notice that the summation of an expression over A or B can be rewritten using its multiplicities.For A this is done like so: Using this equivalence, the function σ, used to calculate the attention weight from x i to x j for A and B, is written as: and Using the condition provided by Equation ( 5), Equation ( 8) can be rewritten as: Substituting back into Equation (7) gives us: The attention weights produced from A can now be expressed in terms of the attention weights produced from B. The j th entry in the attention vector for x i in A is: Thus, the resulting embeddings from the attention mechanism are equivalent.While the embeddings themselves are equivalent, the cardinality for each embedding still differs according to each multisets' multiplicity.The pooling described by Equation (4) fixes this.Let z i be the final embedding for x i .The output vector z from the pooling for A is the sum: Substituting Equation ( 5) again we get, Thus, the output of the PST for both A and B are equivalent.

B Details of MatBench and Jarvis-DFT Datasets B.1 Materials Project
The crystals in the Materials Project [8] are acquired through the use of MatBench [42].MatBench provides standardized training and test sets to ensure that property prediction models can be compared on equal footing.Each training and testing set of the five-fold cross-validation consists of 80% and 20% of the total data respectively.Validation is taken from the training set.In our application we use just 1% of the training data for validation.Due to the relatively simple architecture of the PST, overfitting a target property is more difficult and the typical procedure of selecting the model that performed the most accurately on the validation set is unnecessary.Instead, the PST is trained until the validation error converges.
In Table 7 we list the details of the MatBench dataset for the properties we tested.The information concerning the source of the data and how each of the properties for the crystals therein were calculated as follows.
• Formation Energy -Crystals from this dataset are sourced directly from the Materials Project [8].The property is measured in eV /atom and is calculated using DFT-GGA [50].• Band Gap Energy -Crystals in this dataset are taken from the Materials Project [8].The property is measured in eV and is calculated using DFT-GGA [50].• Bulk/Shear Modulus -Crystals are sourced from the Materials Project.The value of the property is measured in log 10 (GP a).The elastic modulus tensor is derived from the full stress tensor which is calculated using the projector augmented wave method with DFT-GGA functionals [51,50].• Perovskite Formation Energy -Crystals are taken from the work produced by Castelli et al. [52].Property values are measured in eV /unit cell and are calculated with DFT-GGA functionals using the RPBE approximation [53].
• Phonon Peak -Crystals are taken from the Materials Project.Property values are measured in 1/cm and are calculated using DFPT-GGA with PBEsol [54,55].
• Refractive Index -The crystals are taken from the Materials Project.Property values are derived directly from the dielectric tensor [56].Dielectric tensors are calculated using Density Functional Perturbation Theory using VASP with PBE+U exchange-correlation functional [50,57] and Projector Augmented Wave psuedopotentials [51].
• Exfoliation Energy -Crystals are part of the Jarvis-DFT [9] 2D materials dataset, but are originally sourced from the Materials Project and ICSD [58].Property values are measured in meV /atom and are calculated using DFT with optB88 functional [59].Note, that these are different from the crystals mentioned below as part of the Jarvis-DFT dataset for exfoliation energy, which are part of the 3D materials dataset.

B.2 Jarvis-DFT
Much of the Jarvis-DFT database consists of crystals taken from the Materials Project [8] and the ICSD [58].After the initial structures are obtained, they are geometrically optimized using DFT, and then their properties are calculated.Details for each property are as follows: • Formation Energy -The crystals in this dataset are taken from the Materials Project.Their values are measured in eV /atoms and are calculated using OptB88vdW functional [48].
• Band Gap Energy (OPT) -The crystals in this dataset are taken from the Materials Project.The energy values are measured in eV and are calculated using the OptB88vdW functional [48,60].
• Total Energy -The crystals in this dataset are taken from the Materials Project and ICSD.Their values are measured in eV /atom and are calculated using OptB88vdW functional [48].
• Bulk/Shear Modulus -The crystals are sourced from the Materials Project.Property values are measured in GP a and are calculated using the projector augmented wave method using the OptB88vdW functional [61].
• Band Gap Energy (MBJ) -The crystals are taken from the Materials Project.The property values are measured in eV and are computed using the modified Becke-Johnson potential [49,60].
• Spillage -Crystals are taken from the Materials Project.The structures are relaxed using the OptB88vdW functional but Spin-Orbit spillage is calculated using the PBE functional [62].
• SLME -Crystals are taken from the Materials Project.Spectroscopic Limited Maximum Efficiency (SLME) is calculated from the dielectric function and band gap (which are computed using DFT) [9].
• Maximum Piezoelectric Stress Coefficient (e ij ) -Crystals taken from the Materials Project.Value is measured in Cm −2 and is computed by taking the maximum entry from the Piezoelectric stress tensor.
• Maximum Piezoelectric Strain Coefficient (d ij ) -Crystals taken from the Materials Project.Property value is measured in CN −1 and is computed by taking the maximum entry from the Piezoelectric strain tensor.
• Exfoliation Energy -Property values are measured in meV /atom and are calculated using DFT with optB88 functional [9].These structures are different from the crystals mentioned as part of the MatBench dataset for exfoliation energy, which are part of the 2D materials dataset.

C Applications of the Pointwise Distance Distribution
In definition 2.2 we provided a formal definition for the PDD.Here, we will begin by providing an example of its construction.where the first column contains the weights for each row (atom).The second is the distance to the nearest neighbor and the third, is the distance to the second nearest neighbor.The first two rows are identical, as are the final two.Because of this, each of the two groups of rows can be grouped into a single row like so, P DD(S; k) = 0.5 2.481 2.481 0.5 2.881 2.881 The rows are already lexicographically ordered so this is the finalized PDD.If a different unit cell is selected, this collapsing of matrix rows will continue to yield the same result as the proportion of each atom will not change.
The PDDs of two periodic crystals can be compared using the Earth Mover's Distance [63].The weights of each PDD can be considered the distribution we are comparing in the minimum flow problem.We can visualize a set of crystals by projecting these distances onto a two or three-dimensional plane.
Subfigure 4a shows the projection of the pairwise distances of the PDDs of one hundred samples from each of the three crystals in this experiment using Multi-Dimensional Scaling (MDS).This technique attempts to map a set of pairwise distances to specific n-dimensional space by minimizing the difference between the actual pairwise distances and the distances between the points in the projected space.Even with this error, we are able to establish three distinct clusters in two-dimensional space according to their respective crystal type.
Subfigure 4b shows the projection of the pairwise distances between the crystals' PDD onto three-dimensional space for the T2 dataset.The colors of the points in each plot signify the lattice energy for the crystal.The distinction here is not as pronounced; this is to be expected as the sampled crystals are much more similar in their structure and thus the distances between the PDDs are smaller.Nonetheless, there is a discernible trend, and crystals that congregate near each other do share similar lattice energies.

D Prediction of Lattice Energy
Here, the dataset considered contains simulated molecular crystals created by Pulido et al. [65] during the crystal structure prediction using quasi-random sampling [66].This data is subsetted by the underlying molecule, e.g.T2.
During structure prediction, crystals are generated while traversing the potential energy surface and their lattice energy is calculated using DMACRYS [67] to determine their stability.Lattice energy is defined as the energy released during the crystallization of the constituent molecules into a lattice.Prediction of lattice energy is done in three different scenarios.In the first, the model will be applied to a single set of molecular crystals with the same composition using 80%/10%/10% training, validation, and testing splits.Next, the model is applied to multiple sets of crystals each with different underlying molecules using the same splits.The final experiment consists of the application of the model to a set of crystals with an underlying molecule it has not seen in training.The seen data is split 90%/10% for training and validation.To make predictions the PST described in section 2.1 is used with only the PDD as input and no knowledge of the composition.
Invariants are usually used to discern crystals by measuring differences between their structure.Here, the goal is to demonstrate the effectiveness of using an invariant as a representation for a machine learning algorithm.Even when compositional information is not present, the PDD can distinguish crystals with the EMD from the changes in the pairwise distances that occur when the species of atoms are changed.Whether the same distinction can be made in the context of a learning algorithm has previously not been shown.
We make a comparison to another invariant that has been used to predict lattice energy.Average-minimum-distance [33] was used as input for a Gaussian regression model [34].In Table 8 we list the performance on the test set of this AMD model to allow comparison between invariants.
Our model reduces the mean-absolute-error (MAE) rate by 21% compared to the Gaussian Process Regression technique which utilizes AMD [34].The mean-absolute percentage error (MAPE) is also reduced by 1.63%.While we use k = 60 nearest neighbors for constructing the PDD, the AMD model uses k = 500 to achieve its best results.Despite the use of this additional information, the model using the PDD still performs more accurately.In the second row of Table 8 we list the results of the second experiment.The previous task is extended to a dataset of crystals that contains different underlying molecules.These molecules have different compositions.This compositional information is not contained within the PDD (and thus, not in our input).The model will have to discern crystals solely by their structure.While the overall MAE has increased slightly, the percentage error has decreased.The domain on which the lattice energies lie is different for each type of crystal.Using the PDD alone is enough for the algorithm to distinguish the crystal types and predict lattice energy accordingly.
The final experiment uses the data from the P1, P1M, and P2 crystals in the training and validation data.The test set consists of the P2M crystal, which is not part of the either training or validation set.This experiment is the closest to real-world conditions in which new crystals often arise and finding the stable forms is crucial, but information on their lattice energy is unavailable.
When lattice energy is calculated using ab initio calculations, the range of the energies varies from crystal to crystal.When introducing a new type of crystal for our algorithm to make predictions on, this can become a problem as extrapolation to unvisited parts of the lattice energy range can result in high error rates.Fortunately, lattice energies between different types of crystals, are not usually meant to be compared.Instead, they are generated for potential polymorphs of a crystal in an effort to find those with the highest stability for synthesis.We can make use of this fact when applying our model to novel crystals.It is not necessary to predict the correct range of lattice energies; instead, the model needs to be able to predict the lattice energies of the various structures such that their ordering according to their lattice energy is correct.This task could feasibly be turned into a learning-to-rank problem [68], but as a regression task, it allows for a more general approach since the predicted lattice energies can be ordered after the fact.
Each dataset has its lattice energies shifted by the mean lattice energy towards zero.By doing this they each maintain their distribution but now overlap around the origin.The model is trained and validated based on this shifted data.The MAE of the predictions on the test set after they have been shifted back is 7.24 kJ/mol and the MAPE is 3.89%.
While the MAE and MAPE are higher than in previous experiments, the improvement over AMD is more significant.The majority of errors come from underestimating the lattice energy.The datasets tend to grow sparser in these areas where lattice energy is lower as this is where the most stable structures tend to lie.Having a false positive (predicting a higher energy structure to be lower) increases the number of potentially stable structures.False negatives are more impactful as they may result in a structure not being considered entirely due to its seemingly low stability.
The histogram in Figure 5a shows the lattice energies (in kJ/mol) of the three crystals within the dataset.Figure 5b shows the comparison of the predictions of the model against the true property values in the dataset, referred to as ground truth values.Predictions that have lower errors will have their point placed closer to the line colored in blue.The bulk of the points share their error both below and above the true lattice energy as we would expect in a model without bias.There are a few outliers, in particular, a single crystal from the T2 dataset has a predicted lattice energy of just −0.41.Prevention of such errors can be mitigated by using a different loss function than MAE.In particular, mean-squared error (MSE) and Huber loss can hedge against outlying errors by increasing their contribution to the loss function.We choose to not present the results using these loss functions as MAE still provides better results for MAE and MAPE.In Figure 6a we plot the original predictions on P2M after they have been re-scaled back away from the mean lattice energy.In Figure 6b, we re-scale the prediction and the true property values found in the dataset (henceforth, referred to as ground truth values) to between zero and one.By doing this, we can see the ordering of the predictions compared to the true lattice energies.The scatter plot in Figure 6b allows us to compare lattice energies relative to other predictions.
If the error rate produced by the final experiment of section D is inadequate, we can supplement the training with predictions from classical methods.For a novel crystal, ab initio calculations can be used to generate lattice energies for a small subset of structures.These samples can be integrated into the training set to improve the overall results.In order to make this practical, we can only produce lattice energies for a limited portion of structures; we limit this to 10% (or under) of the total structures.9 displays the result of adding this supplemental data to the training set.As expected, the addition of data decreases MAE and MAPE.Even with as few as 367 samples (or 5% of the total dataset), the reduction in error is significant.This process experiences diminishing returns as the amount of the original dataset is used in the training set.Though the MAE and MAPE continue to decrease, it is questionable whether or not spending time generating the instances using classical methods is worth the additional performance gains.

E Additional Experiments E.1 Identifying Crystals with Outlying Property Values
The properties that have been targeted with the PST are all regression tasks.These can also be reframed as a classification task wherein the model attempts to identify if a given material has a "high" (or alternatively, a "low") property value.Due to the ambiguity of this term, a threshold at which a property qualifies as having a "high" value needs to be established.In the following experiments, a value exceeding the 90 th percentile of the test set is considered "high".
MAE is no longer a suitable metric for measuring the model's performance.Instead, four other metrics are used: accuracy, recall, and precision.We define the term true positive (TP) as the correct classification of materials with a "high" property value.False positive (FP) is the incorrect classification of a sample with a "high" property value.
True negative (TN) refers to the correct classification of a crystal that does not have a "high" property value.Finally, False Negative refers to the incorrect classification of material that does not have a "high" property value.Accuracy is determined by the expression (T P + T N )/(T P + T N + F P + F N ).Precision is found by evaluating the expression T P/(T P + F P ) and recall is found using T P/(T P + F N ).The results of this classification task for all properties provided by MatBench are listed in Table 10.The scores listed are the average across all five folds of each dataset.The accuracy for all properties is relatively high, with a lower bound of 94% occurring for the property with the fewest number of samples: exfoliation energy.The trend is similar for recall and precision, with properties containing a smaller sample size having lower scores.Across all properties, the relationship between precision and recall remains balanced with no single property containing an exceedingly high rate of false negatives or false positives.Increasing the threshold at which crystals are classified as positive decreases both recall and precision, but accuracies increase.This can be seen in Table 11 where the threshold is increased to the 95 th percentile.At this level, most properties still perform well with only small reductions in recall and precision.Exfoliation energy is an exception to this.This is the byproduct of having a comparatively small dataset.
If materials with property values above the 99 th percentile are considered positive samples the precision and recall of refractive index drops to 24% and 14% percent respectively.Notably, at this threshold, there are only 7 such crystals which qualify as positive.Datasets with a large number of crystals do not experience as drastic a decrease in efficacy.Band Gap energy, for example, experiences a reduction to 87.1% and 80.4% for precision and recall, respectively.

E.2 Effect of k-nearest neighbors
PDD encoding can be said to be parameterized by two values, the collapse tolerance and the number of k-nearest neighbors.The integer k determines the dimensionality of initial PDD encoding embedding.As k increases, it retains all information from the previous (smaller) values of k.The initial nearest neighbor distances are the most important and the embedding has diminishing returns after this.For any value of k > 1, the PDD is invariant.Thus, if the PDD is different, the crystals are guaranteed to be structurally different.In order for the PDD to be distinct such that if any two crystals are different, their PDDs are different, we need the property of generic completeness.This property is given provided the lattice L and sufficiently large k.An upper bound on this k is when all distances in the last column of the PDD are larger than twice the covering radius of the lattice L of the periodic set.We would expect the performance of the model to increase up until this point.
The upper bound on k can be exceedingly large so we implement a heuristic to find a lower bound that is more computationally efficient.As k increases, the number of rows collapsed in the PDD will either stay the same or decrease.The lower bound on k can be considered an integer large enough that the groups established at the upper bound of k are the same as this lower bound.Each crystal could have a different k for which this requirement is met.
Our encoding method prohibits the use of a dynamic k value, therefore we need a consistent value for k that can be applied to all crystals in the dataset.The results in The error rates listed in Table 12 vary the value of k and report the resulting mean MAE across the five folds.As expected, the lower values of k generally result in higher MAE.Increasing k eventually causes the error rates to stop decreasing.This is also in line with what would be expected as the PDD has enough information to distinguish itself and additional distances are unnecessary.This is not the case only with phonon peak, however, the differences in MAE are relatively small when the deviation of the errors across the folds is considered.

E.3 Effect of collapse tolerance
The collapse tolerance dictates which rows of the PDD will be collapsed.As this parameter increases, the size of the grouped rows will increase.Once rows are grouped, their distances are averaged in the row which represents the group.The change in this averaged row is proportional to the size of the collapse tolerance.In P DD(S; k), as the collapse tolerance approaches infinity, the PDD will decrease in the number of rows until it consists of just a single row with a weight equal to one.In P DD(S; k), the same increase in collapse tolerance will result in a number of rows within the PDD equal to the number of unique elements within the crystal.In both cases, a collapse tolerance that is large enough will result in information loss, eventually increasing errors in predictions.The collapse tolerance is varied and then applied to the Materials Project crystals.The results of this experiment are listed in Table 13.

Property
By increasing the collapse tolerance and reducing the number of rows within the PDD, we can increase the speed of computations.Thus, it is important to choose a tolerance that is maximal, while not sacrificing accuracy.The impact of the collapse tolerance on the size of representation is listed in Table 14.These values are calculated by dividing the number of rows in the PDD by the number of atoms in the unit cell.The number of atoms in the unit cell is used to determine the number of vertices in the crystal graph [12].In this way, the size of our representation can be compared to that of popular graph-based models.Data for crystals typically comes in the form of Crystallographic Information Files (CIF).These files also indicate the amount of potential measurement error for the atomic positions.This is used as a guide and a collapse tolerance of 10 −4 is used in the experiments for the results in Table 1.Sometimes, however, a higher collapse tolerance can act as a regularization technique that is useful on smaller datasets.This effect is seen in the results for refractive index and exfoliation energy, but there is still a balance to be struck.In larger datasets, the performance regression is more noticeable.A collapse tolerance of one is far higher than what would be necessary to cover measurement error in atomic coordinates and would not be advised, even with the potential efficiency gains.Overall, the collapse tolerance does not have a very large impact due to the prevention of rows corresponding to different atoms from being collapsed in the PDD.

Property
Mean

F Implementation Details
The Periodic Set Transformer is implemented using PyTorch [69].There is also a version implemented using Tensorflow [70], however, we have found this version to significantly underperform when compared to the PyTorch version.We believe this to be due to how the output of individual attention head output is handled in their respective implementations of Multi-head Attention.
Data pre-processing is fairly minimal.Each crystal comes in the form of a Pymatgen structure.The structure is converted into a PeriodicSet object.This functionality is provided by the AMD package [33].The PDD of each of the PeriodicSet objects is then calculated with the desired collapse tolerance and k value.Each column is then normalized to between zero and one.This is not necessary for achieving the desired accuracy but it does significantly improve the speed of training by requiring fewer epochs.
With respect to the results in Table 1, training is done on each property with the same hyper-parameters.The hyperparameters that govern PDD encoding remain the same for all properties: a tolerance of 10 −4 and k = 15.Training options including weight decay, epochs, and learning schedule are kept the same as well, except for batch size.Batch size is either 32 or 64 depending on the number of samples: 32 if the number of crystals in the dataset is less than 5000 and 64 if greater.

Figure 2 :
Figure 2: Overview of the architecture of the Periodic Set Transformer.PDD encoding is used to combine the structural information in the PDD with atomic types.The weights of the PDD are incorporated in the attention mechanism and during the pooling of the embeddings to define the multiplicity of the input set.

i
∈ N + in A and B respectively.Each element also carries weight w (a) i ∈ R + and w (b) i ∈ R + for A and B respectively.The application of the Periodic Set Transformer will yield equivalent output if

Figure 3 :
Figure 3: The unit cell of Lutetium-Silicon; Silicon is colored in teal and Lutetium in magenta.

Figure 4 :
Figure 4: Multi-dimensional scaling projection [64] on to R 2 and R 3 for the pairwise distances of crystals between each other.Subfigure (a) projects three types of crystals using distances created by using the Earth Mover's Distance between PDDs for k = 10.Subfigure (b) is created using the MDS projection of pairwise distances between PDDs at k = 100 for one hundred random samples from the T2 crystals colored by lattice energy measured in kJ/mol.
(a) Distribution of ground truth lattice energies.(b) Ground truth lattice energy vs. predictions

Figure 5 :
Figure 5: (a) The distribution of the lattice energies of the T2, S2 and P1 crystals.(b) The predictions of T2, S2, and P1 compared to the ground truth lattice energies in kJ/mol.
(a) Ground Truth vs. Predictions re-scaled according to mean.(b) Ground Truth vs. Predictions normalized between zero and one.

Figure 6 :
Figure 6: (a) Comparison of predictions vs. true (ground truth) values of the test set after the predictions are scaled back by the mean of the lattice energies.(b) The comparison of predictions vs. ground truth after both have been normalized between zero and one.
of the input representation for each dataset in the Materials Project at various collapse tolerances at k = 15 compared to the number of atoms in the unit cell |M |.The size of the input refers to the cardinality of the input set determined by the number of rows in the PDD.The percentage of |M | is the input set's cardinality divided by the number of atoms in the unit cell, expressed as a percentage.

Table 1 :
Five-fold cross-validation prediction MAE and standard deviation of MAE for properties of the crystals in the Materials Project.Bold values indicate the best (lowest) error rate while underlined values indicate the second-best error rate.PST performance is reported using PDD Encoding with a tolerance of 10 −4 and k = 15.encodings can cause ambiguities in the final embeddings, reducing performance.Across all properties the PST outperforms CrabNet, further indicating the usefulness of PDD encoding. the

Table 2 :
[42]le-fold prediction (measured in seconds) and training time (measured in minutes) for the PST (using k = 15 and a collapse tolerance of 10 −4 ) and coGN[22]on Matbench[42].Training and prediction was done using an Nvidia RTX 3090.Time does not include evaluation of the models on the validation sets or data pre-processing times.

Table 3 :
[43]ct of PDD encoding on prediction MAE of the Materials Project crystals.Results are separated by input components where "Composition" uses only the mat2vec atomic embeddings[43]and "PDD" uses only the PDD.Errors in bold indicate the best performance and underlined errors indicate the second-best performance (lower is better ↓).

Table 4 :
Effect of including the PDD weights as defined by Equations (2) and (4) on prediction MAE of the Materials Project crystals.Results for "No weights" use mean pooling and a normal softmax function.Errors in bold indicate the best performance and underlined errors indicate the second-best performance (lower is better ↓).

Table 5 :
[14]iction MAE on the properties of the Jarvis-DFT dataset using the PST and Matformer.Results for Matformer[14]are included for comparison.PST uses PDD encoding with k = 15 and a collapse tolerance of 10 −4 .Bolded values indicate the best performance.The Mean-Absolute-Deviation (MAD) of the test set is included.

Table 6 :
[14]iction (measured in seconds) and training time (measured in minutes) for the PST and Matformer[14]on Jarvis-DFT datasets.Training and prediction was done using an Nvidia RTX 3090.Time does not include evaluation of the models on the validation sets or data pre-processing times.

Table 7 :
Details of the MatBench dataset by property, including the number of samples in each dataset and the meanabsolute-deviation (MAD).

Table 8 :
Results of lattice energy prediction using the PDD with the PST and AMD with Gaussian regression on three different

Table 9 :
Effect of including portions of the P2M data into the training set on prediction accuracy.

Table 10 :
Results of the classification of materials with property values in the 90 th percentile on the MatBench datasets.Threshold refers to the value at which a crystal is determined to be a positive sample if its own property value is higher.The abbreviations in the Table are as follows: (TP) True Positive, (FP) False Positive, (FN)

Table 11 :
Results of the classification of materials with property values in the 95 th percentile on the MatBench datasets.Threshold refers to the value at which a crystal is determined to be a positive sample if its own property value is higher.The abbreviations in the Table are as follows: (TP) True Positive, (FP) False Positive, (FN)

Table 12 :
Table1use k = 15.At this value of k the previously mentioned lower bound is satisfied for 99.1% within the formation energy dataset.Increases in k past this point cause marginal improvements to this coverage that were deemed insufficient when the increased computational cost is considered.Prediction MAE on the Materials Project crystals for various k-nearest neighbors PDD Encoding at a collapse tolerance of exactly zero.Errors in bold indicate the value of k with the best performance and underlined errors indicate the second-best performance (lower is better ↓).

Table 13 :
Prediction MAE on the Materials Project crystals using PDD Encoding at various collapse tolerances with k = 15.Errors in bold indicate the collapse tolerance with the best performance and underlined errors indicate the second-best performance (lower is better ↓).